George Seif

George Seif

library(ggplot2)
library(dplyr)
library(MeanShift)
library(plotly)
library(knitr)
opts_chunk$set(fig.width=9.5)
address <- url("http://www.trutschnig.net/RTR2015.RData")
load(address)
df <- RTR2015[sample(nrow(RTR2015), 10000), ]
plot_ly(df, x = ~longitude, y = ~latitude, size=~rtr_speed_dl)
No trace type specified:
  Based on info supplied, a 'scatter' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
`line.width` does not currently support multiple values.No trace type specified:
  Based on info supplied, a 'scatter' trace seems appropriate.
  Read more about this trace type -> https://plot.ly/r/reference/#scatter
No scatter mode specifed:
  Setting the mode to markers
  Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
`line.width` does not currently support multiple values.
spatial_df_t <- t(df[,c("longitude", "latitude", "rtr_speed_dl")])

Apply standardization and form a set of candidate bandwiths.

spatial_df_t <- spatial_df_t / apply(spatial_df_t, 1, sd)
h.cand <- quantile(dist(t(spatial_df_t)), seq(0.05, 0.40, by=0.05))
system.time(
  bms.clustering <- lapply(h.cand, 
                           function(h){
                             bmsClustering(spatial_df_t, h=h)
                           )
  )

Running blurring mean-shift algorithm...
resulting_df_s <- as.data.frame(t(spatial_df_t)) %>% select(longitude, latitude, rtr_speed_dl) %>% mutate(cluster = as.factor(bms.clustering[[6]]$labels))
(resulting_df <- test_df %>% 
  select(longitude, latitude, rtr_speed_dl) %>% 
  mutate(cluster = as.factor(bms.clustering[[6]]$labels), rtr_speed_dl_s = rtr_speed_dl/sd(rtr_speed_dl)))
g <- list(
  scope = 'europe',
  showland = TRUE,
  # landcolor = toRGB("gray95"),
  # subunitcolor = toRGB("gray85"),
  # countrycolor = toRGB("gray85"),
  countrywidth = 1,
  subunitwidth = 1
)
plot_geo(resulting_df, lon = ~longitude, lat = ~latitude) %>% 
  add_markers(size=~rtr_speed_dl_s
              , color=~cluster
              , text=~paste(paste("Download Speed: ", rtr_speed_dl)
                            , paste("Long: ", longitude)
                            , paste("Lat: ", latitude)
                            , sep = "<br />"
                            )
              , hoverinfo="text"
              ) %>% 
  layout(
    title = 'Mobile download speed in Austria'
    , geo = g
    , mapbox = list(style = "satellite-streets")
  )
`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.`line.width` does not currently support multiple values.
plot( spatial_df_t[1,], spatial_df_t[2,], col=bms.clustering[[5]]$labels, 
xlab="longitude", ylab="latitude", main="Mean shift labels",
cex=spatial_df_t[3,], pch=16 )
points( bms.clustering[[5]]$components[1,], bms.clustering[[5]]$components[2,], col=1:ncol( bms.clustering[[5]]$components ),
pch="+", cex=3 )

p <- plot_ly(spatial_df_t %>% 
               t() %>% 
               as.data.frame()
             , x = ~longitude
             , y = ~latitude
             , z = ~rtr_speed_dl) %>%
  add_markers(color = bms.clustering[[5]]$labels)
p
LS0tDQp0aXRsZTogIk1lYW4gU2hpZnQgQ2x1c3RlcmluZyINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNCiFbR2VvcmdlIFNlaWZdKGh0dHBzOi8vY2RuLWltYWdlcy0xLm1lZGl1bS5jb20vbWF4LzgwMC8xKnZ5ejk0Sl83NmRzVlRvYWE0VkcxWmcuZ2lmKQ0KDQoNCmBgYHtyIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9DQpsaWJyYXJ5KGdncGxvdDIpDQpsaWJyYXJ5KGRwbHlyKQ0KbGlicmFyeShNZWFuU2hpZnQpDQpsaWJyYXJ5KHBsb3RseSkNCmxpYnJhcnkoa25pdHIpDQpvcHRzX2NodW5rJHNldChmaWcud2lkdGg9OS41KQ0KYGBgDQoNCmBgYHtyfQ0KYWRkcmVzcyA8LSB1cmwoImh0dHA6Ly93d3cudHJ1dHNjaG5pZy5uZXQvUlRSMjAxNS5SRGF0YSIpDQpsb2FkKGFkZHJlc3MpDQpkZiA8LSBSVFIyMDE1W3NhbXBsZShucm93KFJUUjIwMTUpLCA1MDApLCBdDQpgYGANCg0KYGBge3J9DQpwbG90X2x5KGRmLCB4ID0gfmxvbmdpdHVkZSwgeSA9IH5sYXRpdHVkZSwgc2l6ZT1+cnRyX3NwZWVkX2RsKQ0KYGBgDQoNCmBgYHtyfQ0Kc3BhdGlhbF9kZl90IDwtIHQoZGZbLGMoImxvbmdpdHVkZSIsICJsYXRpdHVkZSIsICJydHJfc3BlZWRfZGwiKV0pDQpgYGANCg0KQXBwbHkgc3RhbmRhcmRpemF0aW9uIGFuZCBmb3JtIGEgc2V0IG9mIGNhbmRpZGF0ZSBiYW5kd2l0aHMuDQoNCmBgYHtyfQ0Kc3BhdGlhbF9kZl90IDwtIHNwYXRpYWxfZGZfdCAvIGFwcGx5KHNwYXRpYWxfZGZfdCwgMSwgc2QpDQoNCmguY2FuZCA8LSBxdWFudGlsZShkaXN0KHQoc3BhdGlhbF9kZl90KSksIHNlcSgwLjA1LCAwLjQwLCBieT0wLjA1KSkNCmBgYA0KDQpgYGB7cn0NCnN5c3RlbS50aW1lKA0KICBibXMuY2x1c3RlcmluZyA8LSBsYXBwbHkoaC5jYW5kLCANCiAgICAgICAgICAgICAgICAgICAgICAgICAgIGZ1bmN0aW9uKGgpew0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICBibXNDbHVzdGVyaW5nKHNwYXRpYWxfZGZfdCwgaD1oKQ0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICB9DQogICAgICAgICAgICAgICAgICAgICAgICAgICApDQogICkNCmBgYA0KDQpgYGB7cn0NCnJlc3VsdGluZ19kZl9zIDwtIGFzLmRhdGEuZnJhbWUodChzcGF0aWFsX2RmX3QpKSAlPiUgc2VsZWN0KGxvbmdpdHVkZSwgbGF0aXR1ZGUsIHJ0cl9zcGVlZF9kbCkgJT4lIG11dGF0ZShjbHVzdGVyID0gYXMuZmFjdG9yKGJtcy5jbHVzdGVyaW5nW1s2XV0kbGFiZWxzKSkNCihyZXN1bHRpbmdfZGYgPC0gdGVzdF9kZiAlPiUgDQogIHNlbGVjdChsb25naXR1ZGUsIGxhdGl0dWRlLCBydHJfc3BlZWRfZGwpICU+JSANCiAgbXV0YXRlKGNsdXN0ZXIgPSBhcy5mYWN0b3IoYm1zLmNsdXN0ZXJpbmdbWzZdXSRsYWJlbHMpLCBydHJfc3BlZWRfZGxfcyA9IHJ0cl9zcGVlZF9kbC9zZChydHJfc3BlZWRfZGwpKSkNCmBgYA0KDQpgYGB7cn0NCmcgPC0gbGlzdCgNCiAgc2NvcGUgPSAnZXVyb3BlJywNCiAgcHJvamVjdGlvbiA9IGxpc3QodHlwZSA9ICduYXR1cmFsIGVhcnRoJyksDQogIHNob3dsYW5kID0gVFJVRSwNCiAgY291bnRyeXdpZHRoID0gMSwNCiAgc3VidW5pdHdpZHRoID0gMQ0KKQ0KDQpwbG90X2dlbyhyZXN1bHRpbmdfZGYsIGxvbiA9IH5sb25naXR1ZGUsIGxhdCA9IH5sYXRpdHVkZSkgJT4lIA0KICBhZGRfbWFya2VycyhzaXplPX5ydHJfc3BlZWRfZGxfcw0KICAgICAgICAgICAgICAsIGNvbG9yPX5jbHVzdGVyDQogICAgICAgICAgICAgICwgdGV4dD1+cGFzdGUocGFzdGUoIkRvd25sb2FkIFNwZWVkOiAiLCBydHJfc3BlZWRfZGwpDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgLCBwYXN0ZSgiTG9uZzogIiwgbG9uZ2l0dWRlKQ0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICwgcGFzdGUoIkxhdDogIiwgbGF0aXR1ZGUpDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgLCBzZXAgPSAiPGJyIC8+Ig0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICkNCiAgICAgICAgICAgICAgLCBob3ZlcmluZm89InRleHQiDQogICAgICAgICAgICAgICkgJT4lIA0KICBsYXlvdXQoDQogICAgdGl0bGUgPSAnTW9iaWxlIGRvd25sb2FkIHNwZWVkIGluIEF1c3RyaWEnDQogICAgLCBnZW8gPSBnDQogICAgLCBtYXBib3ggPSBsaXN0KHN0eWxlID0gInNhdGVsbGl0ZS1zdHJlZXRzIikNCiAgKQ0KYGBgDQoNCmBgYHtyfQ0KcGxvdCggc3BhdGlhbF9kZl90WzEsXSwgc3BhdGlhbF9kZl90WzIsXSwgY29sPWJtcy5jbHVzdGVyaW5nW1s1XV0kbGFiZWxzLCANCnhsYWI9ImxvbmdpdHVkZSIsIHlsYWI9ImxhdGl0dWRlIiwgbWFpbj0iTWVhbiBzaGlmdCBsYWJlbHMiLA0KY2V4PXNwYXRpYWxfZGZfdFszLF0sIHBjaD0xNiApDQpwb2ludHMoIGJtcy5jbHVzdGVyaW5nW1s1XV0kY29tcG9uZW50c1sxLF0sIGJtcy5jbHVzdGVyaW5nW1s1XV0kY29tcG9uZW50c1syLF0sIGNvbD0xOm5jb2woIGJtcy5jbHVzdGVyaW5nW1s1XV0kY29tcG9uZW50cyApLA0KcGNoPSIrIiwgY2V4PTMgKQ0KYGBgDQoNCmBgYHtyfQ0KcCA8LSBwbG90X2x5KHNwYXRpYWxfZGZfdCAlPiUgDQogICAgICAgICAgICAgICB0KCkgJT4lIA0KICAgICAgICAgICAgICAgYXMuZGF0YS5mcmFtZSgpDQogICAgICAgICAgICAgLCB4ID0gfmxvbmdpdHVkZQ0KICAgICAgICAgICAgICwgeSA9IH5sYXRpdHVkZQ0KICAgICAgICAgICAgICwgeiA9IH5ydHJfc3BlZWRfZGwpICU+JQ0KICBhZGRfbWFya2Vycyhjb2xvciA9IGJtcy5jbHVzdGVyaW5nW1s1XV0kbGFiZWxzKQ0KcA0KYGBgDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQo=